function [P,dP,hP] = logricepdf(x, A, sigma, selGrHess)
% [P] = LOGRICEPDF(X, A, SIGMA)
% Returns the logarithm of the rice probability density of the values X, 
% with amplitude A and a  noise level of sigma. Default: sigma= 1
%
% [P,dP,d2P] = LOGRICEPDF(X, A, SIGMA, selGrHess)
% Additionally compute the first (dP) and second (d2P) derivative of the
% log pdf with respect to the X (selGrHess(1) = true), A (selGrHess(2) = true) 
% and/or, sigma (selGrHess(3) = true). default: selGrHess = [true true true];
%
% Created by Dirk Poot, University of Antwerp, 13-8-2007

if nargin<3
    sigma = 1;
end;
% see http://en.wikipedia.org/wiki/Rice_distribution
rs2 = 1./sigma.^2;
bx = rs2.*x.*A;
besX = besseli(0,bx,1);

%log(AB) = log(A) + log(B). The term +abs(real(bx))) is because the besseli(0,bx,1) 
%function produces a term -abs(real(bx))), which should not occur here.
P = log(x) + log(rs2) +(-.5.*rs2.*(x.^2+A.^2)+ abs(real(bx))) + log(besX);

if nargout>1
    if nargin<4
        selGrHess = true(1,3);
    elseif ~islogical(selGrHess) || numel(selGrHess)~=3
        error('incorrect selGrHess, should be 3 element logical vector.');
    end
    % compute gradient.
    besRat = besseli(1,bx,1)./besX;
    if selGrHess(1)
        dPdx = 1./x + rs2 .* (A .* besRat - x);
    else
        dPdx = [];
    end
    if selGrHess(2)
        dPdA = rs2 .* ( x .* besRat - A);
    else
        dPdA = [];
    end;
    if selGrHess(3)
        dpds = ( (A.^2 + x.^2).*rs2 - 2 * bx .* besRat - 2 ) ./ sigma;
    else
        dpds = [];
    end;
    dP = [dPdx(:) dPdA(:) dpds(:)];
    
    if nargout>2
        % also compute hessian.
        besRat2 = besseli(2,bx,1)./besX;
        dPdAA = rs2 .*( .5* rs2 .* x.^2 .*( 1 - 2* besRat.^2 + besRat2 ) - 1 );
        dPdAs = rs2 .* (2*A - 2 .* x .* besRat + rs2 .*( 2 .* A .* x.^2 .*besRat.^2 - A .* x.^2 .*(1 + besRat2)) ) ./ sigma;
        dPdss = rs2 .* ( 2 + rs2 .* (-3*(x.^2+A.^2) + rs2 .* (x.^2.*A.^2 .*( 2 - 4 *besRat.^2 + 2 .* besRat2)) + 2*A.*x .* (3.*besRat)));
        if isequal(selGrHess,[false true false])
            hP = dPdAA(:);
        elseif isequal(selGrHess,[false true true])
            hP = [dPdAA(:) dPdAs(:) dPdss(:)];
        else
            error('hessian with this selGrHess pattern not yet supported');
        end;
    end;
end;
